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ABSTRACT 

A finite difference scheme for three-dimensional steady 
laminar Incompressible flows Is presented. The Navler-Stokes 
equations are expressed conservatively In terms of velocity 
and pressure Increments (delta form). First order upwind 
differences are used for first order partial derivatives of 
velocity Increments resulting In a diagonally dominant matrix 
system. Central differences are applied to all other terms 
for second order accuracy. The SIMPLE pressure correction 
algorithm Is used to satisfy the continuity equation. Numer- 
ical results are presented for cubic cavity flow problems for 
Reynolds numbers up to 2000 and are In good agreement with 
other numerical results. 

1. INTRODUCTION 

Two basic problems for a finite difference scheme for 
Incompressible Navler-Stokes equations are stability and 
accuracy. Second order accuracy can be achieved by using a 
central difference method, but It Is unstable at high 
Reynolds number. Stability at very high Reynolds numbers can 
be achieved using first order upwind differencing, but It Is 
only first order accurate. Several schemes which combine 
upwind and central differences have been developed [1], but 
these methods are only first order accurate. The second 
order upwind method encounters difficulties with convergence 
and accuracy [2], 

The present approach retains the second order accuracy 
of central differencing while utilizing the stability of 
first order upwind differencing. The incompressible Navler- 
Stokes equations In conservative delta form [3] are formu- 
lated, l.e., momentum equations are written In terms of the 
Increments of primitive variables rather than primitive vari- 
ables themselves. The equations are then decoupled and 
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solved seml-lmpllcltly . By using central differencing for 
the primitive variable terms on the right hand side of the 
equations, the scheme achieves second order accuracy. By 
applying first order upwind differencing for the first order 
partial derivative velocity Increment terms on the left hand 
side of the equations, the method yields a diagonally domi- 
nant matrix system which Is then solved by a line by line 
trldlagonal-matrlx algorithm. To satisfy the continuity 
equation, the SIMPLE algorithm [1,4,5] for pressure correc- 
tion Is used. 

The present finite difference approach Is used to solve 
the cubic cavity flow problem which was studied by numerous 
authors [6-10]. This problem Is Important because It Is a 
three-dimensional steady separated flow with very simple 
geometry and boundary conditions. In this paper, the results 
are compared with those of other numerical methods for 
Reynolds numbers up to 2000. 

2. GOVERNING EQUATIONS 

The Navler-Stokes equations for Incompressible viscous 
flows are: 
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where u, v, and w are velocity components, p Is pressure, 
and Re Is Reynolds number. 


All variables are nondlmenslonal 1 zed by the character- 
istic velocity of the upper surface and the characteristic 
length of the box side. 
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to 


Let uN, v". w", and pN be the quantities corresponding 
u, v, w, and p at time level N. Let 


4u N * u N+1 - u N (5) 

and let 4v N , 4w N , and 4p N be defined similarly. 

Applying the fully Implicit method to Eq. (2) for u N+1 , 
one obtains: 
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After dropping all second order terms and rearranging, 
Eq. (6) becomes 
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The above equation Is decoupled from the others by set- 
ting 4v N , 4 w n , and 4p N equal to zero. With this assump- 
tion, the right hand side of the equation Is exactly the 
steady state Incompressible viscous equation. The accuracy 
of the finite difference method applied to this expression Is 
the accuracy of the scheme. Notice that the usual lineari- 
zation gives the term a/ax (u N 4u N ) Instead of a/ax 
(2u N 4u N ) as above. 
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For completeness, the corresponding expressions for 
Eqs. (3) and (4) are 
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Let u* be an estimate for u N+1 etc., then the 
velocity field (u*, v*, and w*) does not satisfy continuity 
In general. To correct velocity and pressure, the following 
pressure equation, based on the SIMPLE pressure correction 
algorithm [12], Is solved. 
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where ip N Is the pressure correction. Velocity and 
pressure are then corrected as follows: 
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where a p Is an under-relaxation factor [12]. 
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3. NUMERICAL METHOO 


An equally spaced rectangular coordinate system Is used 
throughout this study. For the first order space derivatives 
of the left hand side of Eq. (7), first order upwind 
differencing Is used, e.g., 
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Similar expressions apply to 
and 

Central differencing Is used for all diffusion terms, e.g.. 
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where <p can be any variable. 

All terms on the right hand side of Eq. (7) are discre- 
tized by using central differencing, e.g., 

N N N N 

L_ (AN, . u 1+l .3 ,k u 1+l ,3,k ~ u 1-l,3,k u 1-l,3,k 
ax u ; ' 2 ax 

A line by line trldlagonal-matrlx algorithm Is used. In 
the x-momentum equation, Eq. 7(a), 4u N Is found by sweeping 
Impllclty In the z direction only once. Similarly, 4 v N 
and 4 w n are found by sweeping Impllclty In the x and y 
directions respectively. The Poisson equation for 4p N , 

Eq. ( 8 ), Is swept In all three directions (ADI method). 
Central differences are used for Eqs. ( 8 ) to (12). The 
calculation procedure Is sketched schematically as follows: 
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4. INITIAL AND BOUNDARY CONDITIONS 

For the Initial Iteration, u, v, w, and p are set to 
zero at each grid point, except u = 1 on z = 1 (upper sur- 
face); all values on the boundary remain unchanged, except 
the pressure which Is adjusted by second order extrapolation 
at each Iteration. 

5. RESULTS 

All results presented here use the convergence criteria 
of 0.0001 for maximum $u with an underrelaxatlon factor for 
pressure ap = 0.8. All grid sizes are 30 by 30 by 30 except 
otherwise Indicated. 
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The geometry and the boundary conditions for cubic cav- 
ity flow are shown In Fig. 1. A cubical box with Its upper 
surface moving at constant speed u = 1 Is shown In the 
figure. 

The results for Reynolds number 100 are given In Figs. 2 
to 5. Figure 2 shows the comparison of the velocity profile 
on the vertical center line with those of Cazalbou et al. [8] 
and Goda [6]. The agreement with both results Is excellent. 

In order to show the flow directions for very small 
velocity vectors, all velocity vectors are stretched to con- 
stant magnitude (called velocity field In Ref. [8]). The 
velocity field In the symmetry plane Is Illustrated for 
Re = 100 In Fig. 3. A primary vortex Is clearly shown In the 
figure and a small corner eddy at the bottom corner A, Is 
also visible. Figure 4 shows a larger corner eddy at corner 
A, on the cross section just next to the side wall. 

Figures 5(a) to (d) show secondary flows on cross sec- 
tions from x = 0.63 back to x = 0.47. A pair of large vor- 
tices are observed In reverse flow under the primary vortex. 
The sequence of motion for these vortices are shown from 
Figs. 5(a) to (d). This pair of vortices Is formed from the 
side wall and moves toward the symmetry plane over a very 
short distance. They are mixed and dispersed near the sym- 
metry plane. Another pair of secondary vortices Is formed 
from the side walls close to the moving top surface 
(Fig. 5(c)) and disperses In a very short distance. 

The results for Re = 400, 1000, and 2000 are shown In 
Figs. 6 to 9. The results are compared with other methods In 
Figs. 6(a) to (c). Figure 6(a) shows the excellent agreement 
with the results given by Cazalbou et al. [8] for Re = 400. 
The results given by Goda [6] In the same figure, however, 
are slightly different. For the case of Re = 1000, agreement 
with the results of Cazalbou et al. are very good except for 
only one point near the bottom surface. For Re = 2000, 
excellent agreement with the results given by Cazalbou et al. 
[8] Is shown In Fig. 6(c). 

Velocity vectors In the symmetry plane for Re = 400 are 
shown In Fig. 7(a). The variation of magnitude of the velo- 
city vectors Is shown In this figure. Some of the magnitudes 
of velocity vectors are as small as 0.0002. Figure 7(b) Is 
Identical to Fig. 7(a), except It shows only the direction of 
the velocity vectors. Notice that the eddy at corner A, for 
Re = 400 Is larger than that for Re = 100. In Figs. 7(c) 
and (d), velocity fields In the symmetry plane for Re = 1000 
and 2000 are shown. In these figures, the location of the 
primary vortex center Is closer to the center of the symmetry 
plane. A larger eddy Is clearly seen at the corner B. 

Velocity fields In the vertical plane next to the side wall 
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for Re = 400, 1000, and 2000 are given In Figs. 8(a) to (c). 
They are very different from the case for Re = 100 In Fig. 4. 

A larger pair of upper secondary vortices for Re = 400, 1000, 
and 2000 are observed In Figs. 9(a) to (c). In Fig. 9(c), 
there are three pairs of vortices close to the bottom surface 
Instead of one pair seen for the cases of Re = 400 and 1000. 
It Is similar to the phenomena reported experimentally by 
Koseff et al. [11] for Re = 2000. These vortices only appear 
In a very short distance around x = 0.54. 

The convergence history for Re = 100 Is Illustrated In 
Fig. 10. CPU-time Is only 19.8 sec on a CRAY X/MP computer 
for the 20 by 20 by 20 grid size. CPU-time on a CRAY X/MP 
for different cases are given In Table 1. It shows that the 
present scheme Is fairly fast and efficient for low Reynolds 
numbers. However, when Reynolds numbers are high (greater 
than 1000), the scheme needs some Improvement as shown by 
cpu-tlme In Table 1 . 

6. CONCLUSION 

A diagonally dominant finite difference scheme with 
second order accuracy Is presented. Comparisons of present 
results with other numerical methods are In very good agree- 
ment. CPU-time for the case of Re = 100 with 20 by 20 by 20 
grid points Is only 19.8 sec on a CRAY X/MP computer. 
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TABLE 1. - COMPUTING CPU-TIME 
[CRAY X/MP, error <0.0001. ] 


Reynolds 

number 

Grid size 

AT 

CPU-time, 

sec 

Number of 
iterations 

100 

20 by 20 by 20 

0.1 

19.80 

162 ! 

100 

30 by 30 by 30 

.1 

88.39 

230 

400 

30 by 30 by 30 

.1 

108.71 

286 

1000 

30 by 30 by 30 

.1 

147.97 

393 

2000 

50 by 50 by 50 

.01 

4712.76 

2815 



FIGURE 1. - GEOMETRY AND BOUNDARY CONDITIONS OF CUBIC CAVITY FLOW. FIGURE 2. - COMPARISON OF VELOCITY PROFILES ON THE CENTERLINE 

FOR RE = 100. 
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FIGURE 3. - VELOCITY FIELDS IN THE SYMMETRY PLANE FOR 
RE = 100, AX = 1/30. 
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FIGURE H. - 
THE SIDE 

VELOCITY FIELDS IN THE VERTICAL PLANE NEXT TO 
HALL FOR RE = 100, AX = 1/30. 
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FIGURE 7. - VELOCITY VECTORS IN THE SYMMETRY PLANE. 




(B) RE = 1000/ AX = 1/30. 



FIGURE 8. - VELOCITY FIELDS IN THE 
VERTICAL PLANE NEXT TO THE SIDE 
WALL. 
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